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Abstract 

Background: In this paper, we present a framework for improving the accuracy of fixed-step methods for Monte 
Carlo simulation of discrete stochastic chemical kinetics. Stochasticity is ubiquitous in many areas of cell biology, for 
example in gene regulation, biochemical cascades and cell-cell interaction. However most discrete stochastic 
simulation techniques are slow. We apply Richardson extrapolation to the moments of three fixed-step methods, the 
Euler, midpoint and ^-trapezoidal r-leap methods, to demonstrate the power of stochastic extrapolation. The 
extrapolation framework can increase the order of convergence of any fixed-step discrete stochastic solver and is very 
easy to implement; the only condition for its use is knowledge of the appropriate terms of the global error expansion 
of the solver in terms of its stepsize. In practical terms, a higher-order method with a larger stepsize can achieve the 
same level of accuracy as a lower-order method with a smaller one, potentially reducing the computational time of 
the system. 

Results: By obtaining a global error expansion for a general weak first-order method, we prove that extrapolation can 
increase the weak order of convergence for the moments of the Euler and the midpoint r-leap methods, from one to 
two. This is supported by numerical simulations of several chemical systems of biological importance using the Euler, 
midpoint and ^-trapezoidal r-leap methods. In almost all cases, extrapolation results in an improvement of accuracy. 
As in the case of ordinary and stochastic differential equations, extrapolation can be repeated to obtain even 
higher-order approximations. 

Conclusions: Extrapolation is a general framework for increasing the order of accuracy of any fixed-step stochastic 
solver. This enables the simulation of complicated systems in less time, allowing for more realistic biochemical 
problems to be solved. 

Keywords: Stochastic simulation algorithms, r-leap, High-order methods, Monte Carlo error 



Background modelled by ODEs. For instance, when close to a bifurca- 

Biochemical systems with small numbers of interacting tion regime, ODE approximations cannot reproduce the 
components have increasingly been studied in recent behaviour of the system for some parameter values [6]. 
years, as they are some of the most basic systems in cell Stochastic systems can be modelled using discrete Markov 
biology [1-3]. Stochastic effects can strongly influence the processes. The density of states of a well-stirred stochastic 
dynamics of such systems. Applying deterministic ordi- chemical reaction system at each point in time is given by 
nary differential equation (ODE) models to them, which the chemical master equation (CME) [7,8]. The stochas- 
approximate particle numbers as continuous concentra- tic simulation algorithm (SSA) [9] is an exact method for 
tions, can lead to confusing results [4,5]. In some cases, simulating trajectories of the CME as the system evolves 
even systems with large populations cannot be accurately in time. 

The SSA can be computationally intensive to run for 
alistic problems, and alternative methods such as the r- 
leap have been developed to improve performance [10]. 



w! e l P ^ ^ ^ r, ™, ,™ ■ ,„ realistic problems, and alternative methods such as the r- 
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Instead of simulating each reaction, the r-leap performs 
several reactions at once, thus leaping' along the history 
axis of the system. This means that, unlike the SSA, the 
r-leap is not exact; accuracy is maintained by not allow- 
ing too many reactions to occur per step. The size of each 
timestep, r, determines the number of reactions occurring 
during that step, given by a Poisson random number. 

This gain in speed must be balanced with loss of accu- 
racy: larger steps mean fewer calculations but reduced 
accuracy. Many common r-leap implementations employ 
a variable stepsize, as using the optimal stepsize r at 
each point is crucial for the accuracy of the method [10- 
12]. However a fixed-step implementation can be useful 
in some cases. Although it may be less efficient, it is 
much easier to implement than variable-step equivalents. 
More importantly, the extrapolation framework that we 
describe in this paper requires a fixed-step method. 

The original r-leap as described by Gillespie [10] is 
known as the Euler r-leap, as it can be compared to the 
Euler method for solving ODEs. It has been shown to have 
weak order of convergence one under both the scaling 
r 0 (traditional scaling) [13,14] and V oo (large 
volume scaling) [15], where V is the volume of the system. 
In the same paper, Gillespie also proposed the midpoint 
r-leap method [10], which has higher-order convergence 
in some cases [15,16]. Tian and Burrage [17] proposed 
a variant known as the Binomial r-leap method that 
avoids issues with chemical species becoming negative. 
Only recently has more work been done on constructing 
higher-order stochastic methods. One such method is 
the random-corrected r-leap [18], where at each timestep 
a random correction is added to the Poisson random 
number that determines the number of reactions in that 
step. Given a suitable random correction, the lowest 
order errors on the moments can be cancelled. In this 
way methods with up to weak order two convergence 
for both mean and covariance have been constructed. 
More recently, Anderson and Koyama [19] and Hu et 
al. [20] proposed another weak second-order method, 
the # -trapezoidal r-leap, which is an adaptation of the 
stochastic differential equation (SDE) solver of Anderson 
and Mattingly [21] for the discrete stochastic case. 

In this paper we introduce a framework for improv- 
ing the order of accuracy of existing fixed-step stochastic 
methods by using them in conjunction with Richardson 
extrapolation, a well-known technique for improving the 
order of accuracy of a numerical solver by combining 
sets of simulations with different stepsizes [22] . Extrapo- 
lation was originally developed for ODE solvers but has 
also been successfully applied to SDE methods [23]. Our 
approach has three main advantages: 

(1) It increases the order of accuracy of the methods 
supplied to it. This is desirable for the obvious reason 



that the resulting solutions are more accurate, as well 
as that larger timesteps can be used to reach a certain 
level of accuracy, reducing the computational cost. 
This is discussed further in our Conclusions. 

(2) It can be applied to any fixed-step solver, for instance 
inherently higher-order methods such as the 

0 -trapezoidal r-leap or methods with an extended 
stability region such as stochastic Runge-Kutta 
methods [24]. 

(3) The resulting higher-order solutions can be 
extrapolated again to give solutions with even higher 
order, as there is no (theoretical) limit on the number 
of times a method can be extrapolated (although 
statistical errors can obscure the results if the method 
is too accurate - see Section Monte Carlo error). 

Our extrapolated methods may be useful for researchers 
in biology and biochemistry, as they are easy to imple- 
ment and can accurately and quickly simulate discrete 
stochastic systems that could otherwise be too computa- 
tionally intensive. 

We show how the extrapolation framework can be 
applied to fixed-step stochastic algorithms using the 
examples of the fixed-step Euler r-leap, midpoint r-leap 
[10] and 0 -trapezoidal r-leap [20] methods. The extrapo- 
lation procedure depends heavily on the the existence of 
an appropriate global error expansion for the weak error 
of the numerical method. Once this is known, extrapo- 
lation consists of simple arithmetic. We calculate such 
an expansion for an arbitrary weak first-order method; 
this allows us to use extrapolation in order to obtain 
higher-order solutions. The weak order of all the moments 
of such methods can be improved by extrapolation. To 
reinforce this, we perform a simple error analysis by com- 
paring the equations for the true and numerical mean 
of the Euler r-leap method; we see that its global error 
is order one, and extrapolating it increases the order to 
two for the case of zeroth-order and first-order reactions. 
Using numerical simulations, we demonstrate that this is 
true for two first-order and three higher-order test sys- 
tems with the Euler, midpoint and # -trapezoidal r-leap 
methods. Moreover, the extrapolated methods have con- 
sistently lower errors, and in many cases visibly higher- 
order convergence in the first two moments (the lack of 
convergence in some of the simulations is discussed in 
Section Monte Carlo error). Finally, we demonstrate that 
the extrapolation framework can be used to give even 
higher-order numerical solutions by applying a second 
extrapolation to the Euler r-leap method. 

The rest of this paper is organized as follows. We begin 
with an overview of the SSA and the r-leap methods we 
will use later. We then discuss Richardson extrapolation 
for ODEs and SDEs and introduce the extrapolated dis- 
crete stochastic framework. We give numerical results to 
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support our claims that extrapolation reduces the error of 
fixed-step methods. Finally, we discuss the Monte Carlo 
error and give our conclusions. The derivations of the 
global error expansions for SDEs and discrete stochas- 
tic methods and related material are presented in the 
Appendix. 

Overview of stochastic simulation methods 
SSA 

Gillespie's SSA [9] is a statistically exact method for sim- 
ulating paths from a Markov jump process. The two basic 
assumptions of the SSA are (i) that individual molecules 
are not explicitly tracked, and (ii) there are frequent non- 
reactive collisions. Thus we assume that the system is 
well-mixed and homogeneous. 

The SSA simulates a system of biochemical reactions 
with N species and M reactions, interacting inside a fixed 
volume V at constant temperature. The populations of 
chemical species (as molecular numbers, not concentra- 
tions) at time t are represented as a state vector x = 
X(£) = (Xi, . . . ,Xjv) r . Reactions are represented by a 
stoichiometric matrix vj = (vy, . . . , vxj) T , where ; = 
1, . . . ,M, composed of M individual stoichiometric vec- 
tors. Each stoichiometric vector represents a reaction ; 
occurring and the system changing from state x to x + vj. 
Each reaction occurs in an interval [t,t + r) with relative 
probability aj(x)dt, where aj is the propensity function of 
the ;-th reaction. Propensity functions are given by the 
mass-action kinetics of the reactant chemical species. For 
more detail, the reader is referred to Ref. [9] . The variables 
X, vj and aj(X) fully characterise the system at each point 
in time. 

Algorithm 1. SSA Direct Method 

With the system in state X n at time t n : 

1. Generate two random numbers r\ and r<i from the 
unit-interval uniform distribution U(0, 1). 

2. Find the time until the next reaction r = — In ( — ), 

where a 0 (X n ) = J^ii «/(X„). 

3. Find next reaction ; from 

4. Update t n +\ = t n + r and X w+ i = X n + vj. 

The Direct Method requires two newly-generated ran- 
dom numbers at each timestep. Although there are other 
SSA implementations, such as the Next Reaction Method 
[25] and the Optimised Direct Method [26], which can be 
more economical, in general the SSA is computationally 
costly. 

r-leap method 

The r-leap algorithm leaps along the history axis of the 
SSA by evaluating groups of reactions at once [10]. This 



means significantly fewer calculations, i.e. shorter compu- 
tational time, per simulation, but simulation accuracy is 
compromised: we do not know exactly how many reac- 
tions occurred during each time step, nor can we tell 
more precisely when each reaction occurs than in which 
timestep. The leap condition defines an upper bound for 
the size of each timestep r: it must be so small that the 
propensities do not change significantly for its duration, 
i.e. the change in state from time t to t + r is very small 
[10]. Since r is small, the probability a(x)z that a reac- 
tion occurs during [ t, t + r) is also small, so the number 
of times Kj each reaction fires over one timestep can be 
approximated by V(aj(x)r), a Poisson random variable 
with mean and variance aj(x)z [10]. The Euler r-leap algo- 
rithm is the basic r-leap method, and corresponds to the 
Euler method for solving ODEs or the Euler-Maruyama 
method for solving SDEs. 

Algorithm 2. Euler r-leap method 
With the system in state X n at time t m and a timestep r: 

1. Generate M Poisson random numbers 
kj = V(aj(X n )r),j=l,...,M. 

2. Update t n +\ =t n + x and X n+ i = X n + J2pLi v jkj- 

The Euler r-leap has weak order one [13-15]. Although 
considerable work has been done on improving the mech- 
anism for selecting the timesteps r [10-12] and elimi- 
nating steps that would result in negative populations 
[17,27-29], this does not affect the order of the method, 
limiting its accuracy. Methods with higher order are the 
only way to improve the accuracy beyond a certain point. 
Realising this, Gillespie also proposed a higher-order r- 
leap method, the midpoint r-leap [10]. This is similar to 
the midpoint method for ODEs, where at each step an 
estimate is made of the gradient of X at t n + r /2. X n is then 
incremented using this extra information to give a more 
accurate approximation. 

Algorithm 3. Midpoint r-leap method 
With the system in state X n at time t m and a timestep r: 

1. Calculate X' = X n + \x J^ii vjaj(X n ). 

2. Generate M Poisson random numbers 
kj = P(ajQC)T),j=l,...,M. 

3. Update t n +\ =t n + x and X n+ i = X n + J^fti v jkj- 

Although under the scaling r 0 the midpoint r-leap 
has the same order of accuracy in the mean as the Euler 
r-leap method, under the large volume scaling it has weak 
order two [15,16]. Our numerical simulations also suggest 
that it gives higher-order approximations to the first two 
moments for both linear and non-linear systems (although 
this is not clear from the literature). However the local 
truncation error of its covariance is first-order [16]. 
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0 -trapezoidal r-leap method 

Based on the SDE method of Anderson and Mattingly 
[21], the 0 -trapezoidal r-leap [20] is a weak second-order 
method. It consists of two steps, a predictor step with size 
Or and a corrector step with size (1 — 0)x that aims to 
cancel any errors made in the first step. 
Algorithm 4. # -trapezoidal r-leap method 

For a specified 6 e (0, 1), a\ = 2 (i-6)6 > a2 = 1-0)0 ' 
With the system in state X n at time t m and a timestep r: 

1. Generate M Poisson random numbers 
k'j =V(aj(X n )0r)J= 1,...,M. 

2. Calculate predictor step X' = X n + Y!)Li vfty 

3. Calculate lj = max (a\aj(X!) — a2aj(X n ), 0). 

4. Generate M Poisson random numbers 



kj = V(lj(l-0)T),j=l,...,M. 



5. Update t n +i = t n + r and X w+ i = X! + X!;=i v ; -A/. 

Specifically, the 0 -trapezoidal r-leap method was shown 
to have weak order of convergence two in the moments, 
and a local truncation error of 0(r 3 V -1 ) for the covari- 
ance. r = V - ^, 0 < ft < 1 and in the analysis V — >► oo, 
but in simulations the system volume is kept constant; 
thus it seems that in practice this also results in weak 
second-order convergence in the covariance [20]. 

Methods 

The extrapolation framework 

We start with an ODE solver with stepsize h approximat- 
ing the true solution x( T) by x^ (where T = nh), for which 
we assume the following global error expansion exists: 

x(D-x* = e k (T)h k +e k+1 (T)h k+1 +e k+2 (T)h k+2 +. . . , 

(1) 

where k is the order of the numerical method and the 
ek(T) are constants that only depend on the final integra- 
tion time T. Extrapolating has the effect of cancelling the 
leading error term, resulting in a more accurate approx- 
imation. The existence of such an expansion is key to 
constructing a higher-order approximation, as the appro- 
priate extrapolation coefficients must be used for the 
leading error terms to cancel. For example, in the case of a 
first-order method with stepsize /z, 



x(T) 



e 1 (T)h + e 2 (T)h 2 + (D(h 3 ), 



(2) 



and similarly for stepsize |, 

x(D - 4i 2 = ei(D \ + e 2 (D^ + 0(h\ 

Z 4 

Setting x^ = 2x^ 2 — x^ and using (3) and (2), we obtain 



(3) 



x(T)-i h n = --e 2 (T) + 0(h 3 ), 



(4) 



which implies that x^ is now a second-order 
approximation to x(T). 

h 

We define z(k, q) = x n q to be a series of approximations 
with order k and stepsize h q to the true solution x(T), 
where T = nh\, h q = T/p q and h\ > h 2 > . . . > h q . 
In general, one can use an order k method with step- 
sizes h q -i and h q (in the previous example, h q -\ = h and 
h q — |, i.e. p = 2), to arrive at an order k + 1 estimate to 
x(T), 

x(T)-z(k,h q - 1 ) = 0(h k + 1 ), 



where 



z(k,hq-i) = 



p k z(k, q - 1) - z(k, q) 
p k - 1 " 



This process can be repeated indefinitely. We 
can extrapolate from the initial approximations 
z(l, 1), . . . ,z(l,q) by combining the successive solutions 
in each column of the Romberg table: 

z(l,l) 

z(2,l) 

z(l,2) z(3,l) 

z(2,2) : z(q,q) 

z(l,3) : z(3,q) 

: z(2,q) 
z(l,q) 



For instance, in Eq. (4) we used (with p = 2) x^ = z(l, 1) 
and x^ 2 = z(l, 2) to find x^ = z(2, 1). Repeating with x^ 2 
and x^ 4 , we could extrapolate to find x^ 2 = z(2, 2). Then 
we could extrapolate x^ and x^ 2 to find a third-order 
approximation x^ = z(3, 1), and so on. 

Stochastic methods have two error measures: strong 
(comparing trajectories) and weak (comparing moments). 
In the context of extrapolation, we are interested in the 
global weak error, defined as 



\E{f(x(T)))-W(4))l 



(5) 



where/ : M. N \-> R is a suitable smooth functional, for 
example the first moment of one of the components of x. 
x^ is a numerical approximation to the SDE 



dX t = a(X t , t)dt + b(X u t)dW t , 



(6) 



xM 



where a(x, t) : M iV+1 R n , b(x, t) : h> 
and is a standard M-dimensional Wiener increment. 
Talay and Tubaro [23] derived a similar expansion to Eq. 
(1) for the global error when x^ was calculated using 
the Euler-Maruyama and Milstein schemes (outlined in 



Szekely etal. BMC Systems Biology 201 2, 6:85 
http://www.biomedcentral.eom/1752-0509/6/85 



Page 5 of 19 



Appendix A). By using this expansion and the extrapola- 
tion framework, they were able to derive a second-order 
approximation to K(f(x(T))). The crucial step in obtain- 
ing the global error expansion was to express it as a 
telescopic sum of the local errors. Liu and Li [30] also 
followed a similar procedure to derive a global error 
expansion for numerical methods for SDEs with Poisson 
jumps, thus allowing them to obtain higher-order weak 
approximations. 

Extrapolation for discrete chemical kinetics 

The extrapolation framework can be extended to the dis- 
crete stochastic regime. Since it requires two or more 
sets of approximations with given stepsizes (e.g. h and 
/z/2), it can only be used with a fixed-step method: as 
more complex r-leap methods vary r at each step, it 
is not clear how to extrapolate them. However, this 
has the advantage of making our method very easy to 
program, as there is no complex programming over- 
head, for instance in choosing the timestep for r- 
leap methods. We stress that we mostly use extrap- 
olation to obtain higher-order approximations to the 
moments of the system (or their combinations, such 
as the covariance). In principle, given enough of the 
moments, the full probability distribution at some given 
time could be constructed. This is known as the Ham- 
burger moment problem [31] and in general is a dif- 
ficult problem to solve, as it might admit an infi- 
nite number of solutions. However, in some cases it 
is possible to reconstruct the full distribution from the 
extrapolated moments, as we have a priori knowledge 
about its shape. For instance, when the final distri- 
bution of states is known to be binomial, only the 
mean and variance are necessary for constructing the 
full extrapolated distribution (see Numerical Results, 
System 1). 

In this section we focus on the Euler r-leap method 
(ETL), since this choice simplifies the analysis, but we 
show in Appendix B that any fixed-stepsize method with 
known weak order can be extrapolated. In our numerical 
investigations we show results for the ETL, the mid- 
point r-leap (MPTL) and the 0 -trapezoidal r-leap (TTTL) 
method. Extrapolating the ETL is very similar to extrapo- 
lating an ODE solver. The extrapolated ETL, which we call 
xETL from here on, involves running two sets of S ETL 
simulations for time T = nr. 

Algorithm 5. Extrapolated Euler r-leap method (xETL) 

1. Run S ETL simulations with stepsize r, to get 

sX^, s = 1, . . . , S. 

2. Calculate desired moments 

E s (f(xD) = i TLif(sO- 

3. Repeat steps 1 and 2 using stepsize r/2 to get 



4. Take (m s tf (x^ 2 )) - E s (f (xj))) as the 

extrapolated approximation to the desired moment. 

Algorithm 5 can be easily modified for use with any 
other fixed-step method, by replacing the ETL in Step 1 
with the chosen method. 

It is instructive to use a simple example to see ana- 
lytically the effects of extrapolating the ETL. When 
the propensity functions are linear (i.e. the system 
only contains zeroth-order and first-order reactions), the 
equations for the moments are closed [32,33] and we 
can find explicitly the global error expansion for the first 
moment of our numerical solution (i.e. choose fix) = x). 
The propensity functions can be written as 



[fli(x),02(x),...,0Af(x)] = Cx + d, 



(7) 



where C e R MxN , v eR NxM and d e R , and we define 
W = vC, i.e. W e R NxN . Thus at some timestep m, mz < 
T, the ETL gives (using matrix notation) 



= x T m + vV(r(Cx T m + d)), 



where x T m is an approximation to the true state vector x(t) 
at the m-th timestep (i.e. time t) with stepsize r. Taking 
the expectation of both sides, the ETL evolves the mean as 

E(x£ +1 ) = E(0 + vE(V(r(Cx T m + d))) 

= E(x^) + vE(E(P(r(Cx^ + d))|x^)) (8) 
= (/ + rW)E(x^) + rvd, 

where / is the N x N identity matrix. Note that we can- 
not evaluate the expectation of V(r(Cx T m + d)) directly: 
because x T m is a random variable, we do not know the dis- 
tribution of V{x x m ). Using the law of total expectation we 
can condition on x T m taking a specific value; V(x T m )\x T m 
does have a Poisson distribution. From (8), we see that at 
time T = nx, 

E(x£) = (l+ znW + (^j^W 2 + . . }j n(0) 

+ W~ l (xnW + ^r 2 n 2 W 2 + . . >d. (9) 

The probability density function of the SSA at time t is 
given by the CME [7,8]. The mean fi(t) = E(x(f)) can be 
found from the CME [34]; it evolves as 



dfi(t) 
dt 



Wfi(t) + vd. 



The solution of this is 



li(t) = e wt iL(Q) + e wt I e~ ws vdds. 



Wt I „-Ws^ 



L 



(10) 
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Using a Taylor expansion and basic manipulation, at t = T 
this evaluates to 

fi(T) = + TW + ^T 2 W 2 + . . ^ fi(0) 

+ W' 1 (^TW +^T 2 W 2 + vd. (11) 

Taking (11) minus (9) we see that the global error is 

IL(T) - E(x£) = (\tTW 2 + 0(r 2 )) ^(0) 

+ QrTW + (9(r 2 )^ vd. (12) 

Furthermore the extrapolated error is 

/t(r)-(2E(x^)-E(x^)) = Qr 2 r^ 3 + (9(r 3 )^(0) 
+ Qr 2 TW 2 + (9 3 ^) vd, 

so the leading error term has been cancelled, leaving 
an order two approximation. Such a calculation would 
also apply for the MPTL. The difference is that for lin- 
ear systems the MPTL is second-order convergent with 
respect to the mean [16], and similarly for the TTTL [20]. 
This should be taken into account in order to choose the 
correct extrapolation coefficients. 

The above analysis only applies for the mean of a linear 
system, a very restricted case, but it is useful for demon- 
strating the basic principles of stochastic extrapolation. 
We employ a similar approach to Talay and Tubaro [23] 
and Liu and Li [30] to find a general expression for the 
global error expansion of the moments of a weak first- 
order discrete stochastic method; this is Appendix B. In 
Appendix C we explicitly evaluate this for the particle 
decay system and show that it is equivalent to Eq. [12] 
in this case. Appendix D contains the equations for the 
second moment for the case of linear systems. 

Multimodal systems 

As discussed before, one limitation of our approach is that 
only specific characteristics of the particle distribution can 
be extrapolated, rather than the full distribution. Typi- 
cally we choose these to be the first and second moments, 
as for many systems these are the quantities of interest. 
However, in some cases the moments do not take val- 
ues relevant to the actual dynamics of the system [35,36]. 
This occurs, for instance, in bimodal or multimodal sys- 
tems, which have two or more stable states. Nevertheless, 
our method can be easily generalised to accommodate 
multimodal distributions as follows. 



Algorithm 6. Extrapolated Euler r-leap method (xETL) 
for multimodal systems 

1. Run S ETL simulations with stepsize r, to get 

5 X^,5= 1,...,S. 

2. Plot histograms of the particle populations at time T 
and identify the stable states. 

3. Choose point(s) at which to partition the S 
simulations into p subsets of Si, . . . , S p simulations 
clustered around each stable state. 

4. Calculate desired moments over the subsets of 
simulations, 

Es/(/ 7 (xS)) = 5Eti/GxS),«-=l,.... J p. 

5. Repeat steps 1 and 4 using stepsize r/2 to get 

6. Take (2E S , (f (x^ 2 )) - E Si (f (x£))) , / = 1, . . . ,p as 
the extrapolated approximation to the desired 
moment for each of the p subsets of simulations. 

Algorithm 6 is also simple to code and does not require 
significant extra computational time compared to Algo- 
rithm 5 because the dynamics of the system are found 
from the original simulations that are necessary for the 
extrapolation anyway. The point(s) at which the simula- 
tions are split into subsets can affect the accuracy of the 
results, so must be chosen with some care. In the Numer- 
ical Results (System 5), we apply Algorithm 6 to a bimodal 
system, and investigate the effects of the choice of splitting 
point. 

Results and discussion 

Numerical results 

We simulate some example systems for various stepsizes 
r over time t =[0,T] using three fixed-step numerical 
methods: the Euler r-leap (ETL), midpoint r-leap (MPTL) 
and 0 -trapezoidal r-leap (TTTL) methods (with 0 = 
0.55), and their extrapolated versions, the xETL, xMPTL 
and xTTTL. We plot the absolute weak errors in the mean 
and second moment, i.e. 

|E(x(D) - E s (xJ)|, |E(x(T)x r (D) - E s (x^) r )| 

(13a) 

for the ETL, MPTL and TTTL methods and 

|E(x(r))-2E s (x^ 2 )+Es(x^)|, 

|E(x(Dx T (D) - 2E s (x^ 2 (x^ 2 ) T ) + E s (x*(x*) r )| 

(13b) 

for the extrapolated methods. Here x(T) is the analyti- 
cal solution at time T and Es(/Xx*)) are the moments of 
its approximations given by S simulations of a fixed-step 
method with stepsize r run for n steps. For the linear sys- 
tems, the true solution is calculated analytically; for the 



Szekely etal. BMC Systems Biology 201 2, 6:85 
http://www.biomedcentral.eom/1752-0509/6/85 



Page 7 of 19 



non-linear systems we use the value given by 10 6 or 10 7 
repeats of the SSA (depending on the system). The error 
of a weak order a method with stepsize r is approximately 
Cx a , where C is an unknown constant. To easily see the 
order of accuracy of our results, we plot all the errors 
on log-log plots. Gradients are calculated using a least 
squares fit through the points. The highest level of Monte 
Carlo error, which can be calculated for the linear systems, 
is marked on the appropriate plots as a straight black line. 
Below this level, the absolute error results are, at least in 
part, essentially random (see Section Monte Carlo error). 
We note that in all test systems, the timesteps used were 
all in the useful r -leaping regime: Poisson counts for each 
reaction channel varied between tens to hundreds. 

System 1: Particle decay system 

A simple test system is a particle decay, 



X-> 0, 



: 0.1. 



The initial particle number was X(0) = 10 4 and the sim- 
ulation time was T = 10.4. The units here, and in the 
systems below, are non-dimensional. This system is use- 
ful only as a test problem, but is first-order and easily 
tractable. The analytical mean and second moment are 
-kT 



E(X(T)) = X(0)e~ 
E(X(T) 2 ) =X(0)e 



-kT 



X(0)e 



-2kT 



(X(0)fe- 2kT . 



The average final particle numbers, calculated as above, 
were E(X(10.4)) = 3534.5. We ran 10 8 simulations using 
timesteps r = 0.05, . . . , 0.8. The errors on the mean and 
second moment are shown in Figure 1. In both cases, the 
ETL gives first-order errors and the xETL gives approx- 
imately second-order errors. The MPTL and TTTL also 
converge with second order. The errors of the xMPTL 
and xTTTL are very small, although they do not converge 



with any noticeable order. This is because the values of 
the absolute error for these methods are effectively given 
by their Monte Carlo error, rather than the bias (this is 
a recurring theme in stochastic simulations - see Section 
Monte Carlo error). The maximum level of Monte Carlo 
error was 0.0148 for the mean and 104.8 for the second 
moment. 

In addition, because the final distribution of this system 
is known to be binomial [10], we can construct the dis- 
tribution of the extrapolated solutions from just the mean 
and variance (Figure 2). The dashed lines are the distri- 
butions of ETL simulations with r = 0.05 (blue) and 
r = 0.8 (red), calculated from their histograms, and the 
circles are the full distributions of the xETL using r = 0.05 
(blue) and r = 0.8 (red). The solution of the CME (black 
line) is the true distribution. Both the extrapolated solu- 
tions match the CME solution very well, in both mean and 
overall shape. 

System 2: Chain decay system 

This system consists of five chemical species, each under- 
going a first-order reaction. It forms a closed chain of 
linear reactions, and is intended as a more complicated, 
but still linear, example system. 

Xl % X 2 ^ X 3 % X^h X 5 % X 1} *y = 0.3,/=l,... ,5. 

The initial populations were X(0) = (2500, 1875, 
1875, 1875, 1875) T and simulation time was T = 16. 
Since this system is linear, its expectation and covari- 
ance can be calculated analytically, as shown in Jahnke 
and Huisinga [32]; we used these to calculate the 
true second moment. The average final particle num- 
bers, given by the analytical mean, were E(X(16)) = 
(1971.3, 1996.7, 2025.1, 2020.9, 1986.1) r . We ran 10 8 sim- 
ulations with r = 0.1, . . . , 1.6. Figure 3 shows the absolute 
errors for the mean and second moment of X\. The ETL 
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Figure 1 System 1 absolute errors. Absolute error (1 3) in (a) the mean for the ETL (gradient 1 .0) and xETL (gradient 1 .9), MPTL (gradient 1 .9) and 
xMPTL (gradient 0.5), TTTL (gradient 2.0) and xmL (gradient 1 .5); (b) the second moment for ETL (gradient 1 .0) and xETL (gradient 1 .6), MPTL 
(gradient 1 .9) and xMPTL (gradient 0.5), ^TL (gradient 2.0) and x^TL (gradient 1 .5). The maximum Monte Carlo error levels (black straight lines) 
were 0.0148 (mean) and 104.8 (second moment). Results from 10 8 simulations, each run for T = 10.4. 
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Figure 2 Constructing full extrapolated distributions of System 1 . Distribution of states at time T = 


1 0.4 for ETL with r = 0.05 (blue dashes), 


ETL with r = 0.8 (red dashes), xETL with r = 0.05 (red circles), xETL with 


r = 0.8 


(blue circles). The analytical solution is given by the CME (black 


line). The extrapolated distributions match the CME solution very closely. 







is again approximately weak order one and the xETL 
weak order two. Again, the errors for the MPTL, TTTL, 
xMPTL and xTTTL are very low, although again their 
order of convergence is not quite as high as expected. 
We believe this is due to the unusually high accuracy of 
these methods for closed systems compared to a maxi- 
mum Monte Carlo error level of 0.0132 (mean) and 52.8 
(second moment), which are high relative to the bias of 
these methods. 

System 3: Michaelis-Menten system 

The Michaelis-Menten is a common non-linear test 
system, and represents an enzyme (X2) reacting with a 
substrate (Xi) to make a product (X4). The enzyme and 
substrate form a complex (X3), which can either dissociate 



or undergo a reaction to a product plus the original 
enzyme. It has four chemical species in three reactions: 



Xi+X 2 


ki 


x 3 , 


h 


= 10- 


x 3 


k 2> 


Xl+X 2 , 




= 0.2, 


x 3 


k 3> 


X2 + X4, 


h 


= 0.2. 



Simulation time was T = 16 and the initial popula- 
tions were X(0) = (10 4 ,2 x 10 3 ,2 x 10 4 ,0) r . We used 
10 8 simulations with r = 0.1, . . . , 1.6. There is no ana- 
lytical solution, so in this case we approximated it with 
10 7 SSA simulations. The average final state, given by the 
SSA, wasE(X(16)) = (5927.0, 18716.2, 3283.8, 20789.2) r . 




Figure 3 System 2 absolute errors. Absolute error (1 3) of X] in (a) the mean for ETL (gradient 1 .3) and xETL (gradient 2.4), MPTL (gradient 1 .2) and 
xMPTL (gradient 1 .6), mL (gradient 1 .1 ) and xmL (gradient 1 .9); (b) the second moment for ETL (gradient 1 .3) and xETL (gradient 2.4), MPTL 
(gradient 1 .6) and xMPTL (gradient 1 .7), ^TL (gradient 1 .0) and x^TL (gradient 2.0). The maximum Monte Carlo error levels (black straight lines) 
were 0.01 32 (mean) and 52.8 (second moment). Results from 1 0 8 simulations, each run for T = 1 6. 
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Figure 4 System 3 absolute errors. Absolute error (1 3) of X] in (a) the mean for ETL (gradient 1 .0) and xETL (gradient 2.0), MPTL (gradient 1 .6) and 
xMPTL (gradient 3.3), TTTL (gradient 1 .6) and xmL (gradient 1 .2); (b) the second moment for ETL (gradient 1 .0) and xETL (gradient 2.0), MPTL 
(gradient 1 .6) and xMPTL (gradient 3.3), ^TL (gradient 1 .6) and x^TL (gradient 1 .2). Results from 1 0 8 simulations, each run for T = 1 6. 



The errors in the mean and second moment for X\ are 
shown in Figure 4. The ETL converges with order one 
for 10 8 simulations, and the xETL with order two. The 
MPTL and TTTL have a similar accuracy to the xETL, 
with approximate order two, with the xMPTL and xTTTL 
of approximate order three. 

We investigated the effects of the coefficient of variation 
(CV, standard deviation divided by mean) on this system. 
The CVs of each species, averaged across all r, in System 
3 at T = 16 were CV(X(16)) = (0.01, 0.003, 0.02, 0.004) r . 
In general, a higher CV indicates that the system is more 
noisy. We chose a new set of parameters for System 3 
to give higher CVs: X new (0) = (100, 50, 200, 0) r , /c new = 
(10~ 4 , 0.05, 0.07) r . The CVs using these parameters 
were CV(X new (16)) = (0.05, 0.03, 0.1, 0.07) r , very 
different from the original CVs. However the relative 
errors (absolute error divided by average SSA state) 



at r = 0.1 were very similar: error re i(X(16)) = 
(0.004, 0.0005, 0.003, 0.001) T for the original system and 
error re i(X ne w(16)) = (0.001, 0.0002, 0.0007, 0.002) T with 
the new parameters (note that it is not useful to average 
the errors across all r). This shows that higher CV does 
not necessarily mean higher errors, and the two are indi- 
cators of different characteristics of the system. We have 
focused on the errors as this is the characteristic that we 
want to improve. 

System 4: Two-enzyme mutual inhibition system 

This is a more realistic system involving 8 chemical 
species and 12 reactions [37]. It represents two enzymes, 
Ea and which catalyse the production of compounds A 
and B, respectively. In a classic example of double-negative 
feedback, each product inhibits the activity of the other 
enzyme. For this reason, the system is bistable in A and B: 
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Figure 5 System 4 absolute errors. Absolute error (1 3) of X] in (a) the mean for ETL (gradient 1 .3) and xETL (gradient 2.6), MPTL (gradient 1 .9) and 
xMPTL (gradient 2.3), TTTl (gradient 1 .7) and xmL (gradient 2.7); (b) the second moment for ETL (gradient 1 .3) and xETL (gradient 2.6), MPTL 
(gradient 2.0) and xMPTL (gradient 2.3), ^TL (gradient 1 .7) and x^TL (gradient 2.9). Results from 1 0 7 simulations, each run for T = 3.2. 
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Figure 6 System 5 low peak absolute errors. Absolute error (1 3) of X\ in (a) the mean for ETL (gradient 0.9) and xETL (gradient 0.06), MPTL 
(gradient 0.8) and xMPTL (gradient -0.3), ^L (gradient 0.4) and x^L (gradient 0.0); (b) the second moment for ETL (gradient 0.9) and xETL (gradient 
-0.1 ), MPTL (gradient 0.8) and xMPTL (gradient -0.6), TTTL (gradient 0.3) and xmL (gradient 0.0). Results from 1 0 8 simulations, each run for 7 = 5. 



when there are many particles of A, few particles of B can 
be produced, and vice versa. The reactions are 



We simulated this system for T 
populations of 



3.2 using initial 



E A 


h 


E A +A, 


h 


Eb 


h 


E b +B, 




Ea+B 


h 


EaB, 


h 



E A B + B 
A 

E B +A 
E B A+A 
B 



5 x l(T 4 ,/<4 = 2, 



^ E A B 2 , h = 10" V<6 = 6, 



*8 

hi 

hi 



E B A, 
E B A 2l 



k 7 

fa 



= 5, 



5 X 1(T 4 , kg = 2, 



/cio = 1(T V<ii = 6, 
kn = 5. 



X(0) = (2 x 10 4 ,1.5 x 10 4 , 9500, 9500, 2000, 500, 2000, 500) T , 



where X = (A,B,E A ,E B ,E A B,E A B 2 ,E B A,E B A 2 ) T . We ran 
10 7 simulations of the r-leap using r = 0.005, . . . ,0.08. 
10 6 SSA simulations were used as an approximation 
to the analytical values. The final state of the system 
as given by the SSA mean was E(X(3.2)) = (10420.5, 
4884.4, 3594.7, 1528.0, 4592.7, 3812.6, 3853.8, 6618.2) r . 
The errors for X\ in this system are shown in Figure 5. 
The ETL and xETL gave errors of approximately order 
one and two, respectively. The MPTL and TTTL were 
again approximate order two; the xMPTL had errors very 
similar to those of the xETL, and the xTTTL had very low 
errors of approximate order three. 
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Figure 7 System 5 high peak absolute errors. Absolute error (1 3) of X] in (a) the mean for ETL (gradient 1 .2) and xETL (gradient 0.9), MPTL 
(gradient 1 .6) and xMPTL (gradient 1 .9), ^TL (gradient 1 .0) and x^TL (gradient 0.4); (b) the second moment for ETL (gradient 1 .2) and xETL (gradient 
0.9), MPTL (gradient 1 .3) and xMPTL (gradient 1 .8), TTTL (gradient 1 .2) and xmL (gradient 0.5). Results from 1 0 8 simulations, each run for 7 = 5. 
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System 5: Schldgl system 

The last example we give illustrates how our method 
could work for systems that have multimodal distribu- 
tions. The Schlogl system consists of four reactions [38] 
and is bimodal: 



A + 2X 


h 


3X, 


h 


= 3 x 10 


3X 




A + 2X, 




= io- 4 , 


B 


h 


X, 


h 


= 10" 3 , 


X 


h 


B, 


h 


= 3.5, 



where the populations of species A and B are held constant 
at 10 5 and 2 x 10 5 respectively, so only the numbers of X 
can change. This is a common benchmark system for com- 
putational simulation algorithms. Under certain parame- 
ter configurations, this system has two stable states for X, 
one high and one low. When X(0) is low, the system usu- 
ally settles in the lower equilibrium, and vice versa for high 
X(0). We used X(0) = 250, an intermediate value that 
gives a bimodal distribution. We ran 10 8 simulations until 
T = 5 using r = 0.0125, . . . , 0.2. Because of the bimodal- 
ity, we separated the data into two sets. As the Schlogl 
system is simple enough to be solved using a CME solver 
(see e.g. [39]), we could calculate the density of states of 
X at T = 5. When simulating more complicated systems, 
this can be approximated by a histogram of the distribu- 
tion of species at T (from the initial ETL simulations, for 
instance) to get a reasonable idea of its shape. Average 
final particle number was calculated from the CME to be 
E(X(5)) = 101.3 for the low peak and E(X(5)) = 546.2 
for the high peak. Figures 6 and 7 show the absolute errors 
for this system (low and high peaks, respectively). The 
error levels seem to be similar in both cases. The general 
trend was for the ETL, MPTL and TTTL to converge with 
approximate order one, whereas the extrapolated solu- 
tions did not have a clear order of convergence. It is clear 
that in this case also, the Monte Carlo error is interfering 
with our ability to see a clear order of convergence. 

The gradient of the MPTL mean error seems to 
change from approximately one (Figure 6) to around 1.5 
(Figure 7). It is unlikely that this is due to Monte Carlo 
error, as the error of the MPTL is high enough that this 
should not be an issue. In fact, this is probably due to the 
large volume limit behaviour of the MPTL, discussed after 
Algorithm 3. Because the mean of the high peak is several 
times higher than the mean of the low peak, the system 
is closer to the large volume limit and the weak order of 
the MPTL increases accordingly. Once in the large volume 
limit, the gradient is expected to be two. 

The point at which the data is separated must be cho- 
sen carefully, as it can influence the error results. Figure 8 
shows the distribution of X at T = 5 calculated using a 
CME solver. The three choices of splitting points that we 
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Figure 8 System 5 full distribution, T = 5. The distribution of 
System 5 calculated using a CME solver for T = 5. The three splitting 
points we have compared in this paper are illustrated. 



compare have also been marked. We chose X = 300; this 
is a fairly obvious splitting point in our case. To support 
this, we tested the effects of splitting the data given by the 
ETL at X = 250 and X = 350. We calculated the relative 
change in mean, second moment and variance between 
these two splitting values, averaged over all simulations of 
the ETL and all five timesteps: Table 1 shows that in this 
case the percentage difference is relatively small for the 
mean, becomes larger for the second moment, and is very 
high for the variance. This implies that the choice of split- 
ting point is very important in this case, and should be 
carefully considered. 

Higher extrapolation 

In principle it is possible to use extrapolation in con- 
junction with some fixed-step method such as the 
ETL to obtain increasingly higher-order approximations 
using the appropriate extrapolation coefficients from the 
Romberg table. In practice, the usefulness of repeated 
extrapolations is debatable, as each adds extra computa- 
tional overhead and the higher accuracy can be obscured 
by Monte Carlo fluctuations. However it might be use- 
ful to create up to third or fourth-order methods in 
this way. We tried double-extrapolating the ETL on our 



Table 1 Relative differences in moments for different 
splitting values of Schlogl system 



Moment 


Low peak 


High peak 


Mean 


6.4% 


1 .6% 


Second moment 


21.8% 


2.5% 


Variance 


81.6% 


50.5% 



Relative differences in moments E(f (X)) for data split at X = 250 and X = 350 
for the Schlogl system as percentages of their values when split at X = 300, i.e. 
100* |E(f(X)) S p| it35 o -E(f(X)) S p| it2 5ol/E(f(X)) S p| it3 oo. 
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Figure 9 System 2, double-extrapolated absolute errors. Absolute error of X] for the Euler r-leap (ETL), extrapolated Euler r-leap (xETL) and 
double-extrapolated Euler r-leap (xxETL) algorithm for (a) the mean, gradients 1 .3 (ETL), 2.3 (xETL), 3.7 (xxETL), and (b) the second moment, 
gradients 1 .3 (ETL), 2.4 (xETL), 3.8 (xxETL). Results from 1 0 8 simulations. Double-extrapolation produces consistently higher-order results. 



test systems, with reasonable success. Figure 9 shows the 
ETL, xETL and double -extrapolated Euler r-leap (xxETL) 
errors on species X\ of System 2; in this case we can 
see that the xxETL results are approximately order three 
(or higher). Double-extrapolating gave substantially lower 
errors in almost every system, but it was not always easy 
to determine the order of convergence. A good example 
of this is System 3 in Figure 10. In such systems with rel- 
atively high Monte Carlo error, the approximate solutions 
from the xxETL were obscured by these fluctuations. The 
order of accuracy of the xxETL could be successfully seen 
for most molecular species of System 2, but it was not 
possible for the other systems as this would require a sig- 
nificant increase in the number of simulations, in order to 
further reduce the Monte Carlo error. 

Monte Carlo error 

The weak error we calculate numerically is 

\E(f(x(T)))-E s (f(x r n ))\. (14) 



This error can be separated into two parts 

[E(f (x(D)) - E 00 (/ , (xJ))] + [E 00 (/-(xJ)) - E s (f(xl))] , 

(15) 

where Eqo^x*)) are the theoretical values of the 
moments calculated by an infinite number of simulations 
with stepsize r. The first term is the truncation error of the 
moments from their analytical solutions, i.e. the bias of the 
method, which depends only on the choice of timestep. 
The second term is the Monte Carlo error, which depends 
only on the number of simulations and is given by -^=, 
where C is some constant and S the number of simu- 
lations. The Monte Carlo error can be so large that it 
overwhelms the bias of the underlying numerical method 
completely; in this case all of the numerical results are, in 
effect, incorrect, as they are random fluctuations. 

This formulation is useful when the propensity func- 
tions are linear. In this case, the moment equations are 



a) b) 




Figure 10 System 3, double-extrapolated absolute errors. Absolute error of Xi for ETL, xETL and xxETLfor (a) the mean, gradients 1.0 (ETL), 2.0 
(xETL) and 2.3 (xxETL), and (b) the second moment, gradients 1 .0 (ETL), 2.0 (xETL), 2.1 (xxETL). Results from 1 0 8 simulations. In this case, the xxETL 
errors do not show a clear order of accuracy, but are nonetheless smaller compared to the other methods. 
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closed, so E QO (/'(x^)) can be calculated for the appropri- 
ate numerical method. As an example, consider the mean 
of the ETL: its true value is given by Eq. (10) and the value 
of its numerical approximation can be found by iterating 
Eq. (8). In addition, a similar calculation can be found in 
Appendix D for the second moment. Unfortunately, this is 
not possible for non-linear systems, since in this case the 
equations describing the evolution of the moments are not 
closed any more [40]. 

Such a breakdown of the errors of System 1 with the 
ETL method is shown in Figure 11, using results from 
10 8 simulations. We see that for the case of the ETL and 
xETL, the bias can easily be seen as the Monte Carlo errors 
are relatively low compared to the bias. However, when 
we extrapolate a second time, the bias of the resulting 
estimator is so low that Monte Carlo fluctuations com- 
pletely obscure it, even with 10 8 simulations. This gives a 
poor approximation for the total error. The only way to 
reduce Monte Carlo error is to run more simulations or 
use variance reduction methods. This is a good illustra- 
tion of why it is important to perform a suitable number 
of simulations to get accurate estimates of the moments. 

Variance reduction methods, which aim to decrease the 
Monte Carlo error, are another useful way of reducing 
computational time: because the Monte Carlo error is 
lower, less simulations need to be run for a given accu- 
racy, saving time. This is an important topic in its own 
right and we do not address it in this paper; we refer the 
interested reader to e.g. [41]. It is an active research area: 
recently Anderson and Higham [42] were able to signif- 
icantly reduce the overal computational cost associated 
with the stochastic simulation of chemical kinetics, by 
extending the idea of multi level Monte Carlo for SDEs 
[43,44]. 

Conclusions 

Throughout this paper, we have given reduced compu- 
tational time as a motivation for using the extrapolation 



framework. As this is an important issue, we support here 
our claim that extrapolation speeds up simulations at a 
given error level. Although twice as many simulations 
must be run for a single extrapolation, the loss in com- 
putational time from more simulations is compensated 
for by the significant reduction in error. This is impor- 
tant, as slow runtime is often a limitation of stochastic 
methods. Total computational times and the correspond- 
ing errors in the mean (in brackets) of all three methods 
used in this paper and their extrapolated and double- 
extrapolated versions are shown in Tables 2 and 3 for 
Systems 1 and 4 (10 8 and 10 7 simulations, respectively). 
The estimated time to run the same number of SSA sim- 
ulations is given for comparison. It should be noted that 
all the extrapolated and double-extrapolated r-leap times 
are estimates: the time-consuming part of the extrapo- 
lation method is the two (or more) sets of simulations 
of the original method that must be run; the extrapola- 
tion itself is a fast arithmetic operation. The times for a 
single extrapolation are calculated as runtimeir = h) + 
runtimeir = h/2), and for a double-extrapolation as 
runtimeir = h) + runtimeir = h/2) + runtimeir — h/4). 
The extrapolated methods take several times longer to 
run, but the errors they give are several to hundreds of 
times lower. Most of the exceptions are where the error 
has clearly reached the Monte Carlo level (Tables 2 and 
3, entries marked (-MC)). Besides extrapolation, the other 
obvious way to reduce error is to use a smaller timestep, 
so the real test for the effectiveness of extrapolation is to 
compare the runtimes of cases with similar error values. 
In Tables 2 and 3, we have marked in bold the extrap- 
olated errors that can be directly compared to the base 
method (i.e. read up one row) and take less time to run, 
and similarly for double-extrapolated errors as compared 
to single-extrapolated ones. Although this varies for each 
system and simulation method, the general trend is that 
extrapolation takes less time to give a similar level of 
error, and this pattern was similar for our other example 
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Figure 1 1 Split errors of System 1 using ETL. Split errors for System 1 for (a) the mean and (b) the second moment. The bias can be easily seen 
with the Euler r-leap (ETL) or only a single extrapolation (xETL), but is obscured when we extrapolate a second time (xxETL). Results are from 1 0 8 
simulations. 
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Table 2 Processing times of System 1 
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ETL 


21.3(9.2) 13.3(18.4) 


7.63 (37.1) 


4.08 (74.7) 


2.21 (152.0) 


xETL 


34.3 (0.05) 


21 (0.16) 


11.8 (0.61) 


6.34 (2.6) 


xx ETL 


(-MC) 42.1 (0.015) 


(-MC) 25 (0.016) 


14(0.043) 


MPTL 


21.3(0.024) 13.4(0.070) 


7.65 (0.25) 


4.18(1.0) 


2.24 (4.2) 


xMPTL 


(-MC) 34.7 (0.0088) (-MC) 21.1 (0.0082) 


(-MC) 1 1 .8 (0.0026) 


6.38 (0.041) 


xxMPTL 


(-MC) 42.4 (0.0089) 


(-MC) 25.2 (0.0090) 


(-MC) 14(0.0088) 


mL 


40.8(0.017) 21.1 (0.063) 


1 3.3 (0.26) 


7.65 (1.1) 


4.17(4.2) 




(-MC) 62. 1 (0.00 1 6) (-MC) 34.6 (0.0022) 


(-MC) 20.9 (0.0044) 


11.8 (0.041) 


xxmL 


(-MC) 75.4 (0.0022) 


(-MC) 42.2 (0.0031) 


(-MQ25.1 (0.011) 


Processing times (in thousands of seconds) of 1 0 8 simulations of System 1 using Euler r-leap, midpoint r-leap and 0-trapezoidal r-leap methods and their 
extrapolated and double-extrapolated versions. Absolute errors in the mean from the analytical values are included in brackets. For comparison, implementation of 
1 0 8 simulations of an SSA using the Direct Method would take approximately 1 .2 x 1 0 5 seconds. Entries marked (-MC) are thought to be below the Monte Carlo error 
level; bold entries are faster than the un/single-extrapolated versions with similar error level (i.e. compare with one row above). 


systems. Thus we feel that in general, extrapolation is a 
worthwhile procedure. 

Monte Carlo error is an unavoidable problem when 
using stochastic simulations. The statistical fluctuations 
inherent in stochastic systems can obscure the bias error 
(i.e. order of convergence) of the numerical method if their 
size relative to the bias is large, as the total error is made 
up of these two contributions. A large number S of sim- 
ulations must be run, as the Monte Carlo error scales as 
1/VS. This error varies for each system. Figures 9 and 
10 (and 11) show this clearly: both of the xxETLs have 
total errors of similar size for the same r, but System 2 
has relatively low Monte Carlo error, allowing us to see 
the bias of that system. However, we believe that System 
3 has relatively high Monte Carlo error compared to its 
bias, implying that the xxETL errors we see in the figure 
are all due to statistical fluctuations. It should be noted 


that this seems to happen for all five test problems we use. 
The reason for this is that the extrapolated methods (and 
even the MPTL and TTTL, in some cases) have very high 
accuracy (i.e. low bias error). Since it is only possible to 
run a limited number of simulations, when the bias is very 
small, the total error will be given almost completely by 
the contribution from the Monte Carlo error. 

A contrasting approach to reducing numerical errors 
is the multilevel Monte Carlo method. Originally devel- 
oped for SDEs [43,44], it has recently been extended to 
discrete chemical kinetics [42]. By considering a formu- 
lation of the total error similar to Eq. (15), the multilevel 
Monte-Carlo method aims to reduce it by decreasing the 
Monte Carlo error. Here also many approximate solutions 
are generated with a variety of different timesteps. By 
intelligently combining many coarse-grained simulations 
with few fine-grained ones, it is possible to find a similar 


Table 3 Processing times of System 4 
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Processing times (in thousands of seconds) of 1 0 7 simulations of System 4 using Euler r-leap, midpoint r-leap and ^-trapezoidal r-leap methods and their 
extrapolated and double-extrapolated versions. Absolute errors in the mean of Xy from the SSA values are included in brackets. Running 1 0 7 simulations of an SSA 
using the Direct Method would take approximately 1 .47 x 1 0 7 seconds. Entries marked (-MC) are thought to be below the Monte Carlo error level; bold entries are 
faster than the un/single-extrapolated versions with similar error level (i.e. compare with one row above). 
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level of accuracy to just using fine-grained simulations. In 
contrast, extrapolation uses the same number of coarse 
and fine-scale solutions and gives results which are more 
accurate than the fine-scale solution, by reducing the 
bias instead of the Monte Carlo error. In cases where the 
bias is obscured by statistical errors, using a combination 
of both extrapolation and the multilevel Monte Carlo 
method would be ideal, as it would reduce both sources 
of error. This is an interesting research question and we 
plan to address it in the future. 

In this work, we have extended the extrapolation frame- 
work, which can increase the weak order of accuracy 
of existing numerical methods, to the discrete stochas- 
tic regime. To demonstrate the concept, we have applied 
it to three fixed-step methods, the Euler, midpoint and 
0 -trapezoidal r-leap methods. Thus we have demon- 
strated numerically the effectiveness of extrapolation on a 
range of discrete stochastic numerical methods with dif- 
ferent orders of accuracy for a variety of problems. The 
main requirement to use extrapolation with a numer- 
ical method is the existence of an expression for the 
global error that relates the error to the stepsize of the 
method. Analytically, this is all that must be found to 
show higher weak order convergence of the extrapolated 
method. To extrapolate once, only the leading error term 
need be known; further extrapolation requires knowledge 
of higher terms. We have found the form of the global 
weak error for a general weak order one method; the pro- 
cedure is similar for higher-order methods. This is the real 
power of our approach: it can be applied to any fixed- 
step numerical method. Moreover, further extrapolations 
can raise the order of accuracy of the method indefinitely 
(although beyond a certain point the lower errors will be 
overtaken by Monte Carlo errors). We expect our method 
to be useful for more complex biochemical systems, for 
instance where frequent reactions must be simulated fast 
but accuracy is still important. 

Appendices 

SDE extrapolation 

We summarise below the work of Talay and Tubaro [23] 
on extrapolating SDEs. The global weak error of a stochas- 
tic numerical scheme is 

errorCT, h) = E(f(x(T))) - E(f(x h n )). (16) 

To extrapolate this scheme, we must obtain an expression 
of the form (1). We start with the Kolmogorov backward 
equations, 

u t + Cu =0, 
u(x, T) =/(x), 

where u t = |f , / is a homogeneous function on 
[0,r]x# -> Mandw(x,£) = E(f(x(T))\x(t) = x). 



x are the solutions of the SDE (6) with initial conditions 
x(0) = xo and C is the generator of (6), 

C = a(x, t) • V x + -b(x, t)b T (x, t) : V X V X , 

with A : B = J^ij^ij^ij- Crucially, (16) can then be written 
as 

error(r,/z) = Ew(x 0 ,0) - Eu(x^ T). 

The first term is known; the second can be found by recur- 
sively calculating the error between u when it is evaluated 
from adjacent timesteps: 

n-l 

Eu(x h n , T) = Eu(x 0 , 0)+h 2 Y, Ety(x£, kh) + h 2 K%, 

k=o 

where ^ (x, t) is an error function that is well-behaved and 
specific to each numerical method and K% is a constant of 
unit order. For the case of the Euler-Maruyama and Mil- 
stein methods, the second term is 0(h), so their global 
weak error is [23] 

error(r, h) = -h f E^(x($), t)ds + 0(h 2 ). (17) 
Jo 

Eq. (17) shows that the Euler-Maruyama and Milstein 
methods have global weak order one. It is easy to see that 
extrapolating them leads to solutions with global weak 
order two. 

Discrete stochastic global error expansion 

Here we derive a global error expansion similar to Eqs. (1) 
and (17) for the numerical approximation of moments by 
a weak first-order discrete stochastic method, such as the 
Euler or midpoint r-leap methods. Once the form of the 
error expansion is known, it is clear what extrapolation 
coefficients to use and extrapolation is simple. Similarly to 
the case of SDEs, we start with the Kolmogorov backward 
equations 

d t u + £u = 0 inZ^x[0, T), (18a) 

u(x, t) =f(x) on ill x { r }- ( 18b ) 

Here 

M 

Cu = ^<2y(x) (u(x + Vj) — U(x)) 
7=1 

and 

u(x,t) = E(f(X T )\X t = x), 

with = X(£). Using semigroup notation it is possible to 
formally denote the solution of (18) as 

u(x,t) = e {T ~ t)c f(x), 

which can be useful for quantifying the local error of 
a numerical approximation. By applying a stochastic 
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Taylor expansion [45] to jump processes [16], the one- 
step expansion for the expectation of/ calculated with a 
first-order numerical method should satisfy 



E{f(X h h )\X h 0 =x) = Ahf(x) 



=/(x) + hCf(x) + -Ajf(x), 

j=2 h 



(19) 



where Aj are difference operators associated with the 
numerical method in hand. Furthermore we have for the 
true one-step expansion 
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Figure 12 Local error. Representation of the local error in terms of 
the function u(x, t) = E(f (X^ x )). 



E(f(X h )\X h 0 =x) = e ch f(x) 



,f( x )+h£f(x) + J2- CJ fW< 

j=2 h 



(20) 



Remark 1. An important element in our derivation of a 
global error expansion relates to the boundedness ofu(x, t), 
and its discrete derivative. This boundedess is guaranteed 
[14] when the number of molecules in the chemical system 
is conserved or decreases with time. Proving this in the gen- 
eral case where zeroth-order reactions can add molecules 
to the system is a non-trivial task. One way around this 
problem is to set the propensity functions aj(x) to zero out- 
side a large but finite domain [14]; this is the approach we 
follow here. 

Now if we let e(f, T, h) define the global error in the 
numerical approximation of E(f(Xj)) at time T = nh by a 
first-order numerical method with timestep h we see that 

e(f, T,h) = E(f(X r )|Xo = x) - E(f(X h T )\X h 0 = x) 



= u(x, 0) - E(u(X h T , T)\X h 0 = x) 

N 

= ^E(ii(x£_ 1)A , « ~ DA)) " E(u(X^,tfr)) 



i=l 

N 



= nu(X* h , ih)) - E(u(X&, m, 



i=l 



where we have used the result in Figure 12, and for nota- 
tional simplicity omitted the dependence of the expecta- 
tions on the initial conditions. If we tskegiiy) = u(y, ih) = 
e^-^fiy), then 



E{u{X* h ,ih))-E{u{X\ h1 ih)) = Efe(X^))-Efe(X* )). 



Applying Eqs. (19) and (20) tog if 



e(f, T,h) = J2 E 

i=l 
N 



(^-A 2 )g i (Xt i _ 1)h )+h 3 Rt 



(C 2 -A 2 )e- hC g i . 1 (X h {i _ 1)h ) 



+ fc 3 E(2?f ) 

N U2 



yE [(£ 2 - A 2 )e- hC u(X h (i _ m , (i - l)h)] 



where we have used the fact thatg, = e~ hc gi-\. We thus 
obtain 



N 



e(f,T,h) = ^/z 2 E(^((X^._ 1) ^(i-l)/z)))+/z 3 E( J Rf), 
i=l 

(21) 

where 

1 



y/f e (x, t) = - (£ 2 - A 2 ) u(x, t), 



(22) 



and 



EORf) < C(T), 

from our assumptions on the boundedness of u(x,t). 
Furthemore it is easy to see that 



N 



J2nE\f e ((X h (i _ 1)h) (i-l)h)\ < C(T). 



i=l 
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Using this and results from Talay and Tubaro [23] we 
find that 

r T 

e(f, T,h) = -h / E(f e (X s ,s))ds + G(h 2 ), 
Jo 

where ^e( x > 0 is dependent on the numerical method and 
given by (22). 

An example explicit calculation of the global error 
expansion for a linear system 

In Eq. (12) we calculated the global error on the mean for a 
linear system using the equations for its true and numeri- 
cal solutions. We show how this approach connects to the 
general formulation of the errors given by Eq. (21) by using 
System 1 as an example. We choose this example as it is 
one-dimensional and linear and is simple enough that the 
global error Eq. (12) can be calculated explicitly from the 
general formulation. We define the following notation for 
the discrete difference: 

/»=/(* + v) -/(*), 

where we have assumed that we have only one chemi- 
cal species subject to only one chemical reaction. For this 
particular case, 

C 2 f = a 2 f vv + a 2 a v f vv + aa v f\ 
while for the ETL 

Atf = a 2 f\ 
and thus 

1 2 

\lr e (x, t) = -(a (x)(a(x + v) — a(x))(u(x + 2v, t) 
— 2u(x + v, t) + u(x, t)) 

+ ^a(x)(a(x + v)— a(x))(u(x + v, t) — u(x, t)). 

(23) 

For particle decay, v = — 1 and a{x) = kx. Forg(#) = x, 
the solution to (18) is 

U (x , t ) — X 6 , 

so (23) becomes 



Thus the global error given by the general formulation is 

N 
i=l 

z ,-=i 

= ^E[(l-Kh)S h ]- 1 e-« T E(X») 

2 \ (l- K h)e Kh -l J N 0 

= K -^E(X h 0 ) + O(h 2 ), 

which agrees exactly with what one obtains by calculating 
Eq. (12) for System 1 (i.e. setting d = 0, W = -k). 

Formulae for the second moment of the Euler 
r-leap in the case of linear systems 

We can write the formulae for the numerical and analyt- 
ical evaluation of the second moment of the ETL. This is 
useful for finding the relative contributions of the bias and 
Monte Carlo errors, as well as for explicitly calculating the 
local errors. The ETL evolves the second moment in time 
as 

E(x^ +1 (x^ +1 ) r )=E[(x^ + vP(r(Cx^+d))) 

(x^ + vP(r(Cx^+d))) T ]. (24) 

Using the law of total expectation and writing S T m = 
E(x^(x^) r ), B T m = diag(CE(x^))) and D = diag(d), we 
obtain 

4+i=^ + r^ + r5^ r +r^d r v r +rvd(^) r 

Wwsi t w T Ww^] n d T v T Wvd^i n ) T w T 

+ r 2 vdd T v T + zvB T m v T + rvDv T . 

(25) 

We can now iterate this formula in order to obtain the 
numerical approximation for the second moment of the 
ETL at any timestep. 

The behaviour of the second moment in time as given 
by the CME, S(t) = E(x(£)x r (£)), is [33,40] 

= WS(t) + S(t)W T + fi(t)d T v T + vdfi T (t) 

dt 

+ vB(t)v T + vDv T , 

(26) 

where Bit) = diag(C7*,(£)). By solving Eq. (26) and iter- 
ating Eq. (25), we can use the formulation of Eq. (15) to 
quantify the bias and Monte Carlo errors of the ETL. 
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A similar approach can also be used for the MPTL and 
TTTL. 
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